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Abstract 

The fast computation of the Gauss hypergeometric function 2F1 with all its parameters 
complex is a difficult task. Although the 2-^1 function verifies numerous analytical properties 
involving power series expansions whose implementation is apparently immediate, their use is 
thwarted by instabilities induced by cancellations between very large terms. Furthermore, small 
areas of the complex plane, in the vicinity of 2; = e^*3 ^ are inaccessible using 2F1 power series 
linear transformations. In order to solve these problems, a generalization of R.C Forrey's 
transformation theory has been developed. The latter has been successful in treating the 
2F1 function with real parameters. As in real case transformation theory, the large canceling 
terms occurring in 2F1 analytical formulas are rigorously dealt with, but by way of a new 
method, directly applicable to the complex plane. Taylor series expansions are employed to 
enter complex areas outside the domain of validity of power series analytical formulas. The 
proposed algorithm, however, becomes unstable in general when |a|,|6|,|c| are moderate or 
large. As a physical application, the calculation of the wave functions of the analytical Poschl- 
Teller-Ginocchio potential involving 2F1 evaluations is considered. 

Program Summary 

Title of the programs: hyp_2Fl, PTG_wf 
Catalogue number: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland 

^E-mail address: nmichel@yukawa.kyoto-u.ac.jp 
Phone : 81-75-753-3857 
Fax : 81-75-753-3886 
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Program summary URL: 
Licensing provisions: none 

Computers on which the program has been tested: Intel i686 

Operating systems: Linux, Windows 

Programming language used: C++, Fortran 90 

Memory required to execute with typical data: 

No. of bits in a word: 64 

No. of processors used: 1 

Has the code been vectorized?: No 

No. of bytes in distributed program, including test data, etc.: 
No. of lines in distributed program: 

Nature of physical problem: The Gauss hypergeometric function 2-^1, with all its parameters 
complex, is uniquely calculated in the frame of transformation theory with power series sum- 
mations, thus providing a very fast algorithm. The evaluation of the wave functions of the 
analytical Poschl-Teller-Ginocchio potential is treated as a physical application. 

Method of solution: The Gauss hypergeometric function 2F1 verifies linear transformation for- 
mulas allowing to consider arguments of a small modulus which then can be handled by a 
power series. They, however, give rise to indeterminate or numerically unstable cases, when 
b — a and c—a — b are equal or close to integers. They are properly dealt with through analytical 
manipulations of the Lanczos expression providing the Gamma function. The remaining zones 
of the complex plane uncovered by transformation formulas are dealt with Taylor expansions 
of the 2F1 function around complex points where linear transformations can be employed. The 
Poschl-Teller-Ginocchio potential wave functions are calculated directly with 2-Fi evaluations. 

Restrictions on the complexity of the problem: The algorithm provides full numerical precision 
in almost all cases for |a|, |6|, and |c| of the order of one or smaller, but starts to be less precise 
or unstable when they increase, especially through a, b, and c imaginary parts. While it is 
possible to run the code for moderate or large \a\, \b\, and |c| and obtain satisfactory results 
for some specified values, the code is very likely to be unstable in this regime. 

Typical running time: 20,000 2-P1 function evaluations take an average of one second. 

Unusual features of the program: Two different codes, one for the hypergeometric function and 
one for the Poschl-Teller-Ginocchio potential wave functions, are provided in C++ and Fortran 
90 versions. 

Keywords: hypergeometric, complex analysis, special functions, analytical potentials 
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PACS: 02.30.Fn, 02.30.Gp, 02.60.-x, 02.60.Gf 



Long Write-up 



1 Introduction 

The Gauss hypergeometric function 2-F1 is one of the most important special functions in com- 
plex analysis. Due to its three defining parameters, its structure is extremely rich and contains 
almost all special functions as a particular or limiting case. Important orthogonal polynomials, 
i.e. Chebyshev, Legendre, Gegenbauer, and Jacobi polynomials, can all be expressed with the 
2-Fi function. It is also the case for elementary functions such as power functions, log, arcsin and 
arctan, and more complicated functions such as elliptic integrals. Letting one of its parameters 
go to -l-cxD, one can obtain, as a limiting case, the confluent hypergeometric function of the 
first and second kind, which itself contains, as particular cases, exponential and sine/cosine el- 
ementary functions, as well as Bessel, Hankel, and Coulomb functions. (See Ref. pQ for a recent 
implementation by the first author of the confluent hypergeometric function in the context of 
Coulomb wave functions.) 

Its physical applications are plethora. For example, it is well known that integrals in atomic 
collision descriptions are very often analytically expressed with the 2F1 function [21 El IH El E] • 
It also appears in the Kepler problem, classically perturbed [7| and relativistic [5]. Still in 
the domain of astrophysics, it can serve as a solution in linearized magnetohydrodynamics 
equations [9J, and is found in the study of gravitational waves emission by binary stars, where 
its computation has to be very precise [lOj- The 2-^1 function plays a particularly important 
role in the topic of analytical potentials. Indeed, wave equations seldom admit analytical 
solutions, so that the study of potentials giving rise to wave functions expressible in closed 
form is very much involved. It can provide analytical wave function solutions of the Schrodinger 
equation [HI [121 1131 [El [151 [IB] and of the relativistic Klein-Gordon [17] and Dirac equations 
[T8| [T9l [20] . At the non-relativistic level, the hypergeometric potentials of Natanzon type 
are particularly important due to their six defining parameters, making them very flexible 
[TB I 1 ^ [2 ^ [23]. They contain the Poschl-Teller-Ginocchio (PTG) potential as a particular case 
[2H [25] , which is widely used in nuclear physics due to its fiat shape mimicking the nuclear 
interior [261 [271 [281 [29] . Furthermore, it bears bound, resonant, or scattering states, features 
very interesting for the study of loosely bound or resonant nuclei. As its wave functions contain 
2F1 function components and also because their computation is not trivial, a section of this 
paper and part of the published code are dedicated to this problem as a direct application. 

The numerical computation of the 2-^1 function is a challenging task. Firstly, the power 
series defining 2-^1 has a radius of convergence equal to one, which implies it can be used only 
for z arguments of small modulus. While transformations formulas can minimize the modulus 
of the argument [30j, they cannot be tackled directly in numerical applications because of the 
appearance of canceling diverging terms for b — a and c — a — 6 in the vicinity of integers [31]. 
Moreover, they cannot span all the complex plane, because the points z = e^*^ cannot be 
reached with any of them [321 [33]. Added to that, the case of even moderate |a|, and |c| 
is problematic due to the rapid growth of power series terms with a,b, and c [32]. Analytic 
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continuation of the 2F1 function by use of convergent series in the complex plane has been 
treated in Ref . [34j . where all possible transformations are devised. However, aforementioned 
problems remain in their approach, as numerical problems occur not for integer b — a and 
c — a — b, for which analytic solutions were presented, but for b — a and c — a — b in the 
vicinity of integers. For the case of real arguments, R.C Forrey's transformation theory |31j 
has allowed to remove the divergences induced by near-integer b — a and c — a — 6 at the price, 
however, of very complicated formulas. They rely on Chebyshev polynomial expansions of the 
Gamma function for arguments belonging to [— | : |], real numbers outside this interval being 
handled through the standard Gamma function property T{x + 1) = xT{x). This algorithm 
is, however, clearly impossible to use in the complex plane. The only code calculating the 
2F1 function in the complex plane known by the authors is the one of Numerical Recipes 
[35] . which integrates numerically the hypergeometric equation in the complex plane starting 
from a value of its disk of convergence. While all problems mentioned above are removed 
with this method, it is necessarily slower in practice than direct summation of power series 
and, moreover, demands a careful check of the used complex path of integration, as integrated 
functions should never decrease rapidly on their path to avoid instabilities |35]. Added to that, 
full numerical precision can never be reached with direct integration [35]. In order to have 
both fast and accurate computation of the 2F1 function in the complex plane, the authors 
chose to use R.C. Forrey's transformation theory, but modified so it relies on the mathematical 
structure of the Lanczos numerical approximation of the Gamma function. The Chebyshev 
method problem is then avoided, and formulas are much simpler as well. The missing parts 
of the complex plane, including the two important points z = e^^^ , are precisely handled by 
Taylor expansions of the 2-^1 function, so that the evaluation of 2F1 in the remaining zones of 
the complex plane costs only two normal 2F1 function calls and one quickly converging Taylor 
series evaluation. They can also be handled with non-hypergeometric power series expansions 
[33] , but they are available for non- integer or vanishing b — a only. While this method is very 
efficient for |a|, and |c| of the order of one or smaller, the general case of moderate or large 
|a|, and |c| with z complex remains, however, problematic. 



2 Definition of the Gauss hypergeometric function 

One defines the Gauss hypergeometric function 2F1 through its power series expansion around 
z = 130]: 






(2) 



n=0 




(3) 
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where = 1 ■ x ■ ■ ■ {x + n — 1) is the standard Pochhammer symboL The recursive definition 
of the series term, used in numerical calculations, is also provided in Eqs. (l^|5]) . a and b can be 
arbitrary complex numbers and the series reduces to a polynomial if a G Z~ or 6 G Z~. The 
domain of definition of c depends on a and b: one can have either c G C\Z^, or c G as long 
as a G or 6 G so that c < a or c < b. In this last case, the power series terminates 
before (c)„ can be equal to zero. 

This series always converges for \z\ < 1, conditionally for \z\ = 1, and always diverges for 
\z\ > 1. 2Fi{(i,b,c; z) is defined in the rest of the complex plane by analytic continuation, 
where it is defined everywhere except on its cut on ]1 : +oo[. From now on, we will denote 
2-Fi(a, b, c; z) with its other standard notation F{z), unless necessary. F{z) verifies the following 
differential equation: 

z{l - z)F'\z) + [c-{a + b+ l)z\F\z) - abF{z) = 0. (4) 

In order to have a solution of Eq.(jl]) in the case of undefined F{z), thus with c G Z~, one 

F(z) 

uses the limit of the function — -— for c — > — m, with m G N and T(z) the standard Gamma 

r(c) 

function [30l [36] : 



lim 

c^-m rfc 



(m + 1)! 



(5) 



Eq.(jl]) provides a useful test T of the accuracy of the numerical evaluation of 2-^1 functions 
for z ^ {0, 1}: 



T 



p,,,. ^ \c-{a + b+l)z]F'{z)~abF{z) 
^ ' z(l — z) 



\F(z)\^ + \F'(z)\^ + \F"(z) 



(6) 



F(0) is trivially tested as F{0) = 1, whereas the following Ti test holds for -F(l) as long as 
F(l) and are defined: 



\[c-ia + b + l)]F'{l)-abF{l)l 



|F(l)L+|F'(l)loc + 10- 



-307 



(7) 



3 Basic linear transformation formulas 

In order to be able to calculate F{z) for \z\ > 1, one will follow Ref. [31] by employing linear 
transformation formulas. They have the peculiarity to express F{z) with 2-^1 functions with 
other parameters and arguments a, b, c, and z, so that it is almost always possible to reach 
arguments of modulus small enough so that the 2^1 functions issued from the transformation 
can be evaluated with a power series. The linear transformation formulas used in the code 
consist in [30] : 
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2^1 (a, 6, c; z) 
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(9) 



(10) 



(11) 



+ 



2Fi{c- a,c-b,c- a-b+l;l- z) (12) 



+ 



(13) 




The first four transformations of Eqs. (!8]|9]|10lllip are always numerically stable. Problems, 

however, arise with the transformations of Eqs. fll2lll3p . Indeed, r(2;) is undefined for negative 
integer arguments, so that these two expressions are indeterminate if c — a — 6 for Eg. (1121) or 
6 — a for Eq. fll3p are integers. The case of negative integer c is not important, as these formulas 
are not used in practice for F{z) being polynomial. As 2-Fi(a, 6, c; z) can be perfectly defined 
even if some linear combinations of a, b, and c are integers, this means that the two terms of 
expressions Eqs. (11211131) diverge separately when b — aoic — a — b respectively becomes integer, 
but so that their sums remain finite. This clearly implies that Eqs. (ll2lll3l) cannot be used 
numerically for b — aoic — a — b respectively close or equal to integers due to cancellations 
between the two diverging terms. While it is possible to lift analytically this indeterminacy for 
the integer case, the close vicinity of integers remains problematic [32l [M] . 

In order to solve this problem, we will rewrite Eqs. (11211130 in order to have their indetermi- 
nacy appear in a simple form, so that it can be handled analytically, in the spirit of Ref . [3l] • 
Our method is, however, much simpler than the one of Ref. [31], so that it can be applied 
without problem to the complex case. 

4 Indeterminacy treatment in linear transformation for- 
mulas 

Two kinds of indeterminacy appear while rewriting Eqs. (ll2|13p : a few involving only elemen- 
tary functions, and one the Gamma function. The first case is handled easily using complex 
generalizations of the real C-functions expml(x) = exp(x) — 1 and loglp(a;) = log(l-l-x), precise 
for \x\ ~ 0. They do not exist in Fortran 90, but can be easily coded [2Z]. The first needed 
functions are then: 
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EAz) 



expml(e z) 
e 

rn—l 



n 



for e 7^ , Eo{z) = z, 
z + e + n) I 6no [0;m-l] 



(14) 



, n=0, njtno 



expml 



+ 



m-l ^ N 

^-^ \z + n J 



n=0, nj^riQ 



■ for e ^ 

m— 1 ^ 



sinc(e) 



JJ + n) [0;m-l] + (2^)™ ^ 



. n=0, n^no 

sml vre' 



ri=0, nj^no 



Z + n 



for e 7^ , sinc(O) = 1. 



(15) 



(16) 



where z G C, e G C, m G N, no = — |3?(2;)J (|xj is the closest integer to the real number 
x), and Sno [0;m-i] is equal to one if uq G [0; m — 1] and vanishes otherwise. The expressions 

provided by P^{z) with e 7^ and P^{z) are defined and stable in all cases, as \z + n\ > - 



\/n 7^ 77-0, immediate from previous definitions. Note that P^{z) 



for e ^ 0. 



sine is the standard cardinal sine function, rewritten here for clarity due to its different possible 
definitions. 

The indeterminacy involving the Gamma function is embedded in the following function: 



GJz) 



Gn(z) 



z + e 



T'(z) 



z G 



T{zy ' 

Go{z) = (-1)"+^ n\ , z = -n ,neN. 
The gamma function is very efficiently implemented with the Lanczos approximation 



(17) 



27r I 2; + 7 — 



N 



C0 + J2 



l + i 



where ^{z) > and 7, and the coefficients Cj, i E [0 : N] are chosen so T{z) is precise 
up to numerical accuracy. Although Eq. (fT8l) is valid for 3ft(z) > 0, it is not used close to the 
imaginary axis to avoid numerical instabilities. To handle the rest of the complex plane, the 

TT 

Euler reflection formula r(2;)r(l — z) = — - — - is applied. 

sin(7r2;) 

When e is not too close to zero, i.e. |e|oo > 0.1, a direct calculation of Gt,{z) with Eq. (fT8|) is 
precise. If 2; G Z~ or 2; + e G Z~, and e 7^ 0, there cannot be any numerical cancellation either 

as then ^iz) = or — (2 + e) = 0, and the case of having 2 G Z~ and e = is trivial. In order 



to calculate Ge{z) in the remaining case, prone to numerical error, one introduces the following 
function: 



HJz) 



GJz)T(z + e). 



(19) 



The latter will be evaluated using the fact that the Lanczos approximate expression of the 
Gamma function is elementary. 

From Eq.([T8]), one obtains for < |e|oo < 0.1, and 3?(z) > | or ^{z + e) > i; 



log 



r(^ + 6) 
r(.) 



+ loglp 



V 



loglp 

N 



-2 + 7 



e log(2; + 7 



e - e 



^ {z - 1 + i){z - 1 + i + e) 



N 



z — 1 



i=l 



+ i 



expml ( log 



HJz) 



T(z] 



(20) 



Due to the well-behaved character of expml and loglp for very small values of |e|oo, Eq. (!20i) 
provides a stable implementation of H^. 

The condition < |e|oo < 0.1 is still being verified, but now for 3fJ(z) < | and 3^^(^ + e) < |, 
one has to use the Euler reflection formula along with Eg. (1181) . as well as the intermediate 
function H_^{z + e): 



HJz + e) 



cosivre) + 



r(^) • 

T{z + e)_ 
Tjl-z-e) 
Til-z) 
sin(7re) 



cos(7re) + 



sm 7re 



tan(7r2;) 



tan(7r2;^ 
H.,{z + e) 
1-eHJz + e) 



H^eil -z) + — smc^ 



TT smci^ej 
tan(7r2;) 



(21) 
(22) 



where H^^{1 — z) in Eq. fl2T|) is calculated using Eq. fl20l) and the Euler reflection formula has 
been used to obtain Eq.( l2T]) . Clearly, Eqs. ( 12T|22|) can be precisely evaluated even when |e|oo is 
very small. 

The e = case is immediate using Eq. (!20l) and Eqs. (l21|22l) with e ^ 0: 



Ho{z) 



Hn(z) 



z — 



-2 + 7 



J + log(2; + 7 - -) - 1 



N 

E 

i=l 



z-l + i) 



N 



Co 



1=1 



z-l+i 



HM-z] 



tan(7r^;) 



z < 



(23) 



(24) 



where Hq{1 — z) is calculated using Eq. (l23|) . It has been numerically checked that the loga- 
rithmic derivative of the Lanczos expression (see Eq. lITSi) ). providing Eg. (1251) . was precise up 
to numerical precision to the logarithmic derivative of the Gamma function for 2; G C so that 
^{z) > |. 

Note that to implement sm{TTz) and tan(7r2;) precisely in previous formulas, one has to use 
the integer n = \^{z)\ and calculate instead sin(7r(2; — n)) and tan(7r(2; — n)) due to the finite 
number of digits used for vr. 

Finally, one arrives at the required expression for G^{z) for |e|oo < 0.1, with z G C\Z^ and 
z + ee C\Z-: 

Ge{z) = JJ^^^^ , \z +\n\\oo <\z + t + |m||oo (25) 



T{z + e) 
H.^jz + e) 
T(z) 



k + I'^lloo > 1-2 + e + |"^||oo (26) 



where n = \^{z)\ and m = \^{z + e)J. These two equivalent equalities are used for Gt,{z) in 
order to have the Gamma function ratio occurring in Eqs. fl20ll2ip as small as possible in case z 
or 2; + e are close to a negative integer. 



5 Numerically stable linear transformation formulas 

Eqs.f ll2|[T3|) are now rewritten so they can be implemented without a problem. Let us recall 
that the polynomial case for F{z) is not considered with these two formulas. 



5.1 1 — z transformation formula 

In order to rewrite Eq. f|T2l) . it is necessary for ^{c — a — b) > to be verified. There is, however, 
no loss of generality, as one can apply Eq.([9]) in case 9ft(c — a — h) < to reach the required 
condition. One defines the integer m = |3fJ(c — a — 6)J and the complex number e = c — a — b — m. 
One then obtains for Eq.( |T2l) . after some tedious manipulations. 



F{z) 
A{z) 

Biz) 



smc e 



{A{z)+B{z)) 



Tic) 



m—l 



{a)n{h)n 



e r(c — a)r(c — b) n\ T{1 + n — m — e) 



Tic) 



+00 



n+m 



(27) 
(28) 



e (a + e)„(6 + e)„ ^ [r(a + e)V{b + e)r(n + m + l)r(n + 1 - e) 
(1 - zY {a + e)n+m{b + e) 

n+m 

r(a)r(6)r(n + m + l + e)T{n + 1) 



'l-z 



(29) 



where 2 G C so that |1 — 2| < 1 and e G C*. e = will be considered afterward as a limiting 
case. 



9 



Let us first consider the implementation of tlie finite sum A{z) of Eg. (1281) . If one writes 



A{z) = 0:^(1 — z)^, one obtains: 



n=0 



r(c) 



e r(l — m — e)r(a + m + e)r(6 + m + e) 

Tic) 



-irfm-1)!: 



a. 



r(a + m)T{b + m 
(a + n){b + n) 



0, 



n+l 



q;„ , < n < m — 2. 



(30) 
(31) 



(n + 1)(1 — m — e + n 

+ 00 

In order to implement B{z), one writes B{z) = /3„(1 — 2;)", and using Eq. (l29ll one obtains 

{a)m{b)m 



n=0 



X 





r(l - e)r(a + m + e)T{b + m + e)T{m + 1) r(a)r(6)r(m + 1 + e) 
r(c)(l-2)" . 



) I >- 1 00 
1 



> 0.1 



r(a + m + e)T{b + m + e) \T{m + 1 
1 



GJl] 



+ GJm + r 



G^{a + m) ^ Ge{b + m) 



r(m + 1 + e) \V{b + m + e) V{a + m) 
E,{\og{l-z)) 



r(a + m)T{b + m)T{m + 1 + e) 



7o 



T{c){a)m{b), 



zy 



r(a + m + e)T{b + m + e)T{m + l)r(l - e 



X r(c)(a)^(b)„(l - 2)™ , |eU < 0.1, (32) 

(33) 



X 



{a + m + n + e){b + m + n + e) 
(m + n + 1 + e) (n + 1) 

{a + m + n) — {b + m + n) — e + 

In 



Pn + 



{a + m + n){b + m + n) 



(m + n + 1) 
a + m + n + e){b + m + n + e) 



n + l 



In+l 



(n + m + l + e)(n + l — e 
{a + m + n){b + m + n) 



n > 0, 



(n + m + l)(n + 1 — e) 



7n , > 



(34) 
(35) 



where 7^ is an auxiliary term needed for the calculation of 

One can see that the expressions defining A{z) and B{z) through their respective series 
terms a„ and /3„ are advantageous for numerical computation. Firstly, Gamma and elementary 
function evaluations only occur in the first term of the series, while all other are rational 
functions of their predecessors. This implies a fast evaluation of the series. Secondly, the 
complication induced by the treatment of small |e|oo is also present only at the level of the first 
term of A[z) and B{z), so that it is not prohibitive to be in the small |e|oo regime compared to 
the normal regime. 
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5.2 - transformation formula 

z 

The method used for the transformation in Eq. (fT3|) is very similar to the one of Sec. (15.ip . In 
this case, the relation 9^(6 — a) > must be verified. Here again, there is no loss of generality, as 
Eq-® can be applied in the case 3fJ(6 — a) < 0. We will use the same notations as in Sec. fl5.ip . 
except for m = |3fJ(6 — a) J and e = b — a — m. One then has for Eq. (fT3i) : 



F{z) 
A{z) 

B{z) 



-1)" 



smc(ej 

r(c) 



iAiz) + B{z)) 



(a)„(l - c + a)r 



£ T{b)T{c - a) ^ r(n + l)r(l + n - m - e) 



(36) 
(37) 



r(c) 



+ 00 



(a)„+m(l - c + a) 



n+m 



e ^ Lr(a + e)r(c - a)V{n + m + l)V{n + 1 - e)(a + e)„(6 + e) 



(-2;) ' (a + e)„+m(l - c + a + e) 



n+m 



(a + e)„ r(a)r(c - a - e)r(n + m + 1 + e)r(n + 1) 
where 2; G C so that Ul > 1 and e G C*. 



z 



— (n+m) 



(38) 
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Their series terms read: 



r(c) 



(3o 



7o 

Pn+l 



+ 



X 



7n+l 



e r(l — m — e)r(a + m + e)r(c — a 

'-^)'"""-^" r(a + ,')r(c-a 

(a + n) (1 — c + a + n) 



6^0 

e = 0, 



(n + 1)(1 — m — e + 



ttn , < n < m — 2, 



r(c)z- 



(a)m(l - c + a)r 



r(c - a)r(a + m + e)r(l - e)T{m + 1) 



(1 - c + a + e), 



leL > 0.1 



r(a)r(c-a-e)r(m + l + e) 

(1 _ c + a + e)„G_,(l) - (1 - c + a)i(l - e) 



{a)mT{c)z~ 



+ (1-c + a + e) 



r(c — a)T{a + m + e)r(m + 1) 
G,{m + 1) Ge{a + m) 



r(c — a)r(a + m + e) r(c — a)r(m + 1 + e) 
G^e{c — a) — i?_e(log(— 2;))p(c — a — e) 



r(m + 1 + e)r(a + m) 
(a)^(l - c + a)^T{c)z-^' 



, |e|oo < 0.1, 



r(a + m + e)r(c - a)r(l - e)T{m + 1) ' 
(a + m + n + e)(l — c + a + m + n + e) 
(m + n + 1 + e)(n + 1) 
(a + m + ra)(l — c + a + m + n) 



(m + n + 1) 
(1 — c + a + m + n) — e + 

7!^ 

(n + m + 1 + e) + 1 — e) 
{a + m + n){l — c + a + m + n) 



— [a + m + n) 



[a + m + n + e)(l — c + a + m + n + e) 



n + l 



(n + m + l)(n + 1 — e) 



, n > 0, 

7„ , n > 0. 



(39) 
(40) 



(41) 
(42) 



(43) 
(44) 



Obviously, this scheme has the same numerical advantages as the 1 — z transformation 
formula of Sec. fl5.ip . 



6 Power series convergence test 

In the formulas of Sec. ([3]), the power series defined in Eq.([T]) is used with various sets of 
parameters. Due to the presence of Pochhammer symbols, its general term modulus may 
increase as n goes to +oo, before finally going to zero. The problem to tackle is the possible 
appearance of false convergence. If one truncates the power series when its general term modulus 
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is smaller than a given precision, one may obtain a wrong or inaccurate value if it eventually 
appears to increase afterward. In Ref. [31], the rest of the power series was majored by an 
analytical series in order to evaluate the needed number of terms for convergence. In order to 
avoid a false convergence, it is not necessary, but sufficient, to know whether the modulus of 
the general term of Eq. ([1]) is either increasing or decreasing to zero for n large enough. 

For that, one defines a polynomial P{x) built upon the ratio of the modulus of two consec- 
utive general terms (see Eq.([n])): 

P{x) = \z {a + x){b + x)f - \{c + x){x + l)f (45) 

which is clearly a polynomial of degree four in x when a, b, c, and z are fixed. P{n) has the 
property to be positive when the general term of Eq. ([T]) increases, while it is negative when it 
decreases. The knowledge of the sign of P{x) would demand the calculation of its four roots, 
which would be too lengthy even with analytical formulas. It is sufficient, however, to consider 
its first and second derivatives, the latter being a polynomial of degree two. Indeed, with A the 
discriminant of P"{x), A < implies that P{x) has only one maximum while A > means that 
it has two local maxima and one local minimum. If A < and P'{xo) < for a real xq, one will 
have P'{x) < Vx > xq. Thus, false convergence becomes impossible to occur once an integer 
no verifying P'{no) < is reached. If A > 0, one considers the largest real root Xc of P"{x), 
and He the smallest integer larger or equal to Xc- The curvature of P{x) becoming negative for 
X > Xc, P{x) has at most one maximum in this zone. Hence, a similar property is obtained as 
for A < 0, that if P'(xo) < for xq > Xc, P'{x) < OWx > Xq, so that false convergence cannot 
appear with integers n > ric- Consequently, during the power series evaluation, it is sufficient 
to calculate P'{n) as long as P'{n) > for n > (A < 0) or for n > ric {A > 0), and, when 
P'in) < 0, to check if its general term at integer n is smaller than numerical precision. There, 
the power series can be truncated safely. As this condition arrives in general very quickly, this 
method does not add expensive additional test calculations which would slow the power series 
evaluation. If Eqs. fll2fl3l) are used, this test is considered for the two power series present in 
these transformations. 



7 Numerical algorithm 

We consider first the non-polynomial case, implying that c G C\Z~. 
7.1 General case 

The cut in [1 : +oo[ is handled replacing z >l\)y z + iO~ (limit from below the real axis), so 
that z becomes z — 2l0~^°^ in the code. One begins with applying the linear transformation of 
Eq.([8]) if 3?(6 — a) < and the one of Eq.([9]) if 3fJ(c — a — b) < in order to be able to to use 
Eqs. (ll2fl3p . The latter transformation also has the advantage to accelerate Eq.([T]) convergence, 

as its term behaves like ^-a-h ^"^^ ^ ~^ +oo, easy to demonstrate using Stirling's approximation 
for Pochhammer symbols and factorial. 

If |1 — 2;| < 10~^, the transformation of Eq. (fT2!) with the method of Sec. (15. II) is applied, as 
it handles the vicinity of the singular point z = 1 properly. 
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Then, the 



and 



runs from 0.5 to 0.9. 



z- 1 
R 



moduh are compared respectively to a given radius i?, which 



0.9 may appear seemingly large, but it was acknowledged that 
power series summations are still fast and precise in this regime if |a|, \h\ and |c| are small or 
moderate. Moreover, although the number of terms necessary to converge may be smaller using 
Eqs.f ll2|[T^ . this advantage is hindered by the larger number of multiplications and divisions 
per series term as well as the presence of several Gamma function evaluations. 
z 

happens to be smaller than R, Eq.([T]) or Eqs. ffT|TU]) are processed, respec- 



If 



z\ or 



1 



tively. Other transformations are considered thereafter if \z\ > 0.9 and 



1 



> 0.9. 



The 



11 — z\, and 



1 



1- z 



moduli are compared respectively to R, which runs 



again from 0.5 to 0.9. If one of them is found to be smaller than R, the corresponding transfor- 
mation, i.e. the one using Eq. f|T3l) . Eqs. fll0|13l) . Eq. f|T2|) . or Eqs.f fT0l[T2l) . respectively, is applied. 
Note that the conditions 9^(6 — a) > and 3ft(c — a — 6) > are always verified even when the 
transformation of Eq. flTUl) is utilized. Eq. flT^ . however, is considered only if |a|oo < 5, \b\oo < 5, 
and |c|oo < 5. Indeed, the transformation of Eq. (fT2!) has been found to become quickly unsta- 
ble when |a|, \b\, and |c| increase. Transformations involving Eq. fllOl) are also processed only if 
|c — 6|oo < 5 in order to avoid instabilities. The number of terms needed to reach convergence 
with the power series is determined with the method explained in Sec. OH]). 

If none of the aforementioned transformations succeed, it has been chosen to determine 
F{z) through a Taylor series expansion. Indeed, the missing zones of our linear transformation 
method, located in the vicinity of z = e^*3 ^ verify 0.9 < \z\ < 1.1, so that for small or moderate 
|a|, and |c|, they are close enough to accessible regions to be tackled with such an expansion 
without inducing instabilities. For this, one defines the center of the Taylor expansion disk 



To—, where Tq = 0.9 if \z\ < 1 and Tq = 1.1 if \z\ > 1. One then has: 



with: 



+00 

n=0 



go = F{zo) , qi = F'{zo) = —2^1(0 + 1, 6 + 1, c + 1; Zo) 

c 



(46) 



qn+2 - 

n> 



Zo{l- Zo){n + 2) 



/ 1N / 7 .NX {a + n){b + n) 
[n{2zo - 1 - c + a + 6 + l)zo) qn+i + -qn 



(47) 



where go and gi are calculated with Eq.([T]) if ro = 0.9, and with Eq. (l36|) if tq = 1.1. No 
expression in denominators may induce false convergence in the three-term recurrence relation 
of Eq. (H7I) . contrary to Eq.([3]), where c + n can become very small and then make the general 
term of the series increase abruptly in modulus. On the contrary, here g^ behaves smoothly. 



so that testing convergence with |g„(-2 — -2o)"|oo + \qn+i{z — Zq) 
instances of very small g„, is sufficient. As \z 



n+ll 



to avoid eventual single 



-2^0 1 ^ 0.1, the Taylor series of Eq. fH^ converges 
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very quickly, with 10-20 terms typically. Hence, virtually all computational time will be spent 
calculating F{zq) and F'{zq), for which zq lies close to the circle of convergence. The zones in 
the vicinity of 2; = e^* 3 thus demand twice as much time as in the general case. This technique 
has been compared to the scheme devised in Ref . [33] for b — a not close to integer. Indeed, the 
latter demands that 6 — a is not a strictly positive integer, and hence not a complex number 
very close to an integer as well, z is taken to be equal to re* 3, with r = 0.99 or r = 1.01, 
so that ^(2^0) and -^'(-^^o) are implemented with respectively Eq.([T]) and Eq. flSB]) . Results are 
shown in Tabs. ffTP]) . If r = 0.99, both the method of Ref. [55] and ours have comparable times 
of calculations, respectively ~ 8 seconds and ~ 10 seconds to calculate Tab.([T]), because Eq.([T]) 
is very quickly summed and no Gamma functions are evaluated therein. However, the method 
of Ref. [33] is faster for r = 1.01 by a factor of 4-5, as it also takes ~ 8 seconds to calculate 
Tab.([l]) with it, whereas it takes ~ 35 seconds with Eqs. (14611471) . Nevertheless, it is clear from 
Tabs. (|T1l2l) that our method is more robust when \a\, \b\ and |c| are not small. Hence, due to its 
ability to treat b — a equal or close to integers properly and due to its robustness, our method 
is to be preferred to that of Ref. [33] . even if it is slower for \z\ > 1. 

7.2 Polynomial case 

Z 

If F{z) is a polynomial, either a G Z~ or 6 G Z~. If a G Z~, one uses Eq.([l]) if |2:| < 

z 1 

and Eqs. (IT][TU]) if not. As both related power series terminate, it is possible to use them with 
\z\ > 1. Indeed, the power series argument modulus cannot exceed two, so that the evaluation 
of the polynomials remains stable. It was checked that for a, b, and c parameters of a given 
order of magnitude, this method provides the same or better precision as the more involved 
non-polynomial case (see Sec. (|TTl) ). If b E Z~ , the procedure is the same, using Eq. (|TT]) instead 



8 Application: wave functions of the PTG potential 

Potentials sustaining analytical wave function solutions of the Schrodinger equation are not 
often. Among them, the class of Natanzon potentials [16j is widely studied as they depend on 
six parameters, so that they can accommodate various shapes. The PTG potential belongs to 
a subclass of the Natanzon potentials and depends on four parameters. Its interest is twofold: 
firstly, parameters can be selected so that it develops a flat bottom, a feature appreciated 
in nuclear theory as it models the interior of nuclei (see Refs. [261 [27] [28 ] for examples of 
applications). Secondly, it possesses a quasi-centrifugal barrier, behaving properly for r — >• 0, 
although vanishing exponentially quickly for r — > +oo, as well as the possibility to bear an 
effective mass. 

8.1 PTG Schrodinger equation 

The Schrodinger equation allowing for the PTG potential reads [21]: 



of Eq.([TU]). 



r 



d_J_d_ ^ £(£+ 1) 
dr /i(r) dr r'^^{r) 



) 



+ VpTcir) u{r) = e u{r) 



(48) 



2mo 
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with rriQ the particle free mass, //(r) its dimensionless effective mass (the full effective mass is 
mo A*(r)), i its orbital angular momentum, Vptg{^) the PTG potential, detailed in this section, 
u{r) its wave function, and e its energy, r has the dimension of a length and niQC^ (c is the 
speed of light), as well as e, of an energy. 

VpTG depends on four parameters, conventionally denoted as A, s, u, and a. A controls the 
shape of the potential so that it is very diffuse for small A, becoming the usual Poschl- Teller 
potential for A = 1, while for larger values of A, a flat bottom develops, s is a scaling parameter 
as VpTG is proportional to its square, and p determines the depth of the potential. Finally, 
a wields the strength of the effective mass and varies between and 1, (1 — a) mo being the 
effective mass at r = 0, while the free effective mass mo is recovered for r — > +oo. 

/i(r) and Vptg{'>^) are expressed through the variable y depending on r, defined by way of 
an implicit equation |24j : 

A^s r = arctanh(?/) + y/A^ - 1 arctan(VA2 - 1 y) , A > 1 (49) 
= arctanh(2/) - Vl - A2 arctanh(Vl - A2 y) , A < 1. (50) 

A numerical technique to quickly and precisely solve this equation will be delineated afterward. 

VpTG^r) is separated in three parts: V^(r), directly proportional to a, Vf(r) its ^-dependent 
part, and Vc{r), its principal central part. One can now outline the expressions for fi{r) and 

vpTGir) m- 

.2\ 



VAr) 



VJr) 



VpTG^r) 



[l-a + [a(4 - 3A') - 3(2 - A')] y^ 
" {l-y')[l + {K'-l)y'], 

l-2/2)(l + (A2-l)y2) 



iii + l) 



-A^iy(u + l) - 



{A'-l)m-a) + 2ay') y'] 



r > 0, 



A' 



(2 - (7 - A V - 5(A2 - 1)?/^) 



(V^ir) + Veir) + V;(r)) 



(51) 

(52) 

(53) 
(54) 
(55) 



2mo/i(r) 

The expression of Ve{r) in Eq. fl53|) presents numerical cancellations for r ^ 0, as y ~ sr 
therein so that both terms of Eq. fl53l) diverge but their sum remains finite. Divergence is 
suppressed if one employs Eqs. fl49ll50p to rewrite Eq. fl53l) so that Eq. fl53|) can be evaluated via 
a power series: 



Vf(r) 



SA{y) 



E- 



y 



-A^sr \ sr 
;i - A2)"+2 



l + -)5A(y) + A^-l-(l + (A^-l)2/' 



r > 0, 



n=0 



2n + 3 
A2-2 



-y 



211 



(56) 
(57) 

(58) 



Eg. (1561) will be used when r > and parameters of arctanh and arctan in Eqs. fl49f5UI) are 
smaller than 0.01, so that Eq. fl57l) will converge quickly. 1^(0) is analytical and provided by 
Eq. flSHl) . With this method, virtually full accuracy is obtained for Vf(r) Vr > 0. 
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8.2 y(r) function calculation 

The parameter y{r) defined in Eqs. fl49|50p has to be handled numerically, as it will be expressed 
analytically only for the Poschl- Teller potential with A = 1, where it is equal to tanh(A^sr). 
Two problems appear while solving this equation. The first one is related to the rapid variations 
of arctanh(y) for y ~ 1, which occurs when r increases. The second problem comes from the 
saturation of y to the value y = 1 for r — »• +oo. When ?/ ~ 1 then 1 — y"^ quickly becomes 
inaccurate by a direct computation. A precise knowledge of 1 — j/^ is indeed demanded to 
calculate properly /i(r) and Vptg(^) (see Eqs.( !5T|55!) ). 

Although y{r) is not analytical, it is possible to locate it analytically inside the [0 : 1[ 



interval. Indeed, y belongs to the interval [yd : yd, where yd and ye read: 

yd = max(tanh(A2sr - VA"^ - 1 arctan(VA2 - 1)),0) , A > 1 (59) 

= tanh(A2sr) , A^ < 1. (60) 

ye = tanh(A2sr) , A^ > 1 (61) 

= tanh(A2sr + Vl - arctanh(Vl - A2)) , A < 1. (62) 

Moreover, y[r) verifies simple expansions for y ^ and y ~ 1: 

y{r) = sr + 0{y^) , r ^ (63) 

= + O ((1 - yf) , r ^ +00 , A > 1 (64) 

= ye + 0{il- yy) , r -> +00 , A < 1. (65) 



As the behavior of y at the extremities of [0 : 1[ is known, it is possible to have a quickly 
converging Newton method procedure to find y{r), as one can there find a very good starting 
point, and also because the derivative, according to y, of Eqs. fl49|50l) is analytical. One uses 



the following starting point ystart for the Newton method: 

Vstart = ydiiyd>^ and A > 1 (66) 

= ?/e if ?/e > ^ and A < 1 (67) 

= sr (0.99 if sr > 0.99) for other cases . (68) 



It is complemented by the bisection of [yd : ye] for the rare cases where y falls out of the interval 
[0 : 1[ during the Newton procedure. 

This procedure, however, cannot be used if y ~ 1 (in the code, the condition yd > 0.99 is 
employed), as numerical inaccuracies would occur. Hence, another method has to be employed. 
For this, one uses an iterative fixed-point scheme on x = arctanh(2/): 

Xn = A^sr - VA2 - 1 arctan(VA2 - ly„) , A > 1 (69) 
= A^sr + Vl - A2 arctanh(Vl - A^y^) , A < 1, (70) 
yn+i = tanh(xn) , n > 0, (71) 

and 2/0 = 1 as a starting point. As 2/ ~ 1, Xn converges quickly to x and y„ to y ior n +oo. 

4g-2a; 

This representation allows a precise determination of 1 — y^, as 1 — = ^ — -2xy ' expression 
not subject to numerical inaccuracy. 
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8.3 PTG wave functions 

8.3.1 Bound, antibound, and resonant states momenta and energies 

The PTG potential sustains bound, antibound, and resonant states, i.e. S-matrix poles, whose 
energies are known analytically [23]. The nature of the state depends on the value of the integer 
N = 2n + i + 1, where n is its principal quantum number, and of the following boundary values: 



Ba = AWtttt^ — T { ly +-]--, A> (72) 

^' A2 1 - a -2 V 2/ 2 ' A/ 1 . V ; 



+oo,A<^/-^, (73) 
1 — a 

= [u] M\N, (74) 

B^ = u-1 ,ueN. (75) 

Bound states occur when < Ba and < B^, antibound states for A^ < Ba and N > By 
[N > if e N), and resonant states for A^ > Ba- Their linear momenta and energy read: 

1 — a 

e = ^ (77) 

where the ± sign is "+" for bound and antibound states whereas it is "-" for resonant or 
complex virtual states, and A reads: 



A = A2(z. + iy(l-a)-[(l-a)A2-l] 



-4' 



2 



(78) 



N = u is suppressed for G N in the bound/antibound case as it provides the value k = 0, 
for which it is easy to show that the PTG wave function vanishes identically. 

8.3.2 Wave functions analytical expressions 

The general expression of PTG wave functions is very similar for S-matrix poles and scattering 
states [23]: 

ip{r) = Af x{r)X+{r)F{r) (79) 
= Af xir){A+X+{r)F+{r)+A-X-{r)F-{r)) (80) 

where A/" is a normalization constant depending on the nature of the state and used functions 
read: 



X^r) = (x+)t , X-{r) = (x+)-f , x(r) = J^^^==^(^-)'*', (81) 

F{r) = ,F,(^iy-,u+,i+^;x~y (82) 
F+(r) = 2F1 {u-,iy+, 1 + P; x+) , F-{r) = 2F1 (//", 1 - P; x+) . (83) 
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Their arguments and employed parameters are: 




^^=1 U^+o -/3MA2(l-a)-l), (85) 



- r(/.+)r(^-) - r(.+)r(.-) ' ^''^ 

= .g^. 

l-y2 + A2^2'-^ l-2y2 + AV 

where the r dependence is directly visible in the and x~ arguments via the y parameter of 
Eqs.f l49ll50ll . From previous expressions, one can check that Eq.( |80ll is obtained by operating 
on Eq. fl79|) with the linear transformation formula of Eq. f|T2|) . The alternative provided by 
Eqs.f l79|80|) for the evaluation of ip{r) is necessary as Eq. flTOj) and Eq. flHOj) are respectively 
unstable for large and small r. Eq. flHOl) is also the usual expression of the wave function in 
terms of outgoing ("+" terms) and incoming part ("— " terms), so that A~ = for S-matrix 
pole states. Occurrences of integer b — aorc — a — b values in 2Fi{a,b,c; z) expressions of 
Eqs.( l82]l83l) can be verified to be possible for Lp{r) scattering (see Eqs. fl79f8Qp ) and for a discrete 
set of linear momenta A; G C. 

Scattering states are normalized according to the Dirac delta normalization. The M nor- 
malization reads: 



_ / 2K\sP{l + \ + B + 2n)r(£ + | + ^ + n)T{l +\ + n) 

^' (£ + I + M'(l - a) + 2n)T{n + l)r(^ + n + l)T{l + |)2 ""'^ 



^ — , ' / \' scattering states . (89) 

27rr(/3)r(-/5)r(£ + |)2 v & ; \ ) 

Asymptotic expressions for y9(r) for r — > and r +oo are obtained from Eqs. fl63ll64|[65l) 
and Eqs. fl79]l80|) : 

ip{r) ~ M v/A(l-a)(Asr)^+^ = C7o/+^ , r 0, (90) 
(^(r) ~ M (A+e^'=(''-"i) + A-e-*^("~"i)) = C+e^'^'" + Q- e-'^"" , r ^ +oo (91) 



with: 



ri = ^ i^K^ - 1 arctan(VA2 - 1) - log , A > 1 (92) 



(^Vl - A2 arctanh(Vl - A^) + log , A < 1. 



(93) 



The code also provides the first and second derivatives of the PTG wave functions, express- 
ible with functions as well. 
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9 The program hyp_2Fl 

The code hyp_2Fl is written in two independent versions, standard C++ and Fortran 90. It 
uses only standard libraries and is thus portable on many machines. 

The C++ code is separated into three different files: complex Junctions. hyp_2Fl.cpp and 
hyp_2Fl_example.cpp in order to ensure compatibility with the code of Ref.p^. Indeed, to use 
this code and Ref. [T] simultaneously, one has to use the present complexJunctions.H instead 
of the file provided in Ref.[l]. The Fortran 90 code for 2F1 is situated in hyp-2Fl.f90 and in 
hyp-2Fl_example.f90. 

complex Junctions. H contains elementary complex functions which are not in the standard 
library, and routines calculating constants specific to the Coulomb wave functions for compat- 
ibility with Ref.lpLj. The only differences with Ref. [I] are situated in log_Gamma, calculated 
with full numerical precision instead of 10"^'' precision, and in the added function GammaJnv, 

which calculate precisely —{z). Only the functions of the complex Junctions. H file needed for 

the 2F1 algorithm appear in hyp-2Fl.f90, to which elementary functions which are present in 
C++ standard libraries but absent in Fortran 90, are added. 

In hyp_2Fl.cpp and hyp_2Fl.f90, routines calculating the 2F1 function consist of: 

• Gamma_ratio_diff_small_eps: calculate H^{z) for < |e|oo < 0.1 (see Eq. fll9p ). 

• GammaJnv-dijf-cps: calculate Ge{z) (see Eq.l ITTll ). 

• A_sum_init: calculate — — of the ao term of A(z) (see Eqs.( l30|39l) ). 

e r(l — m — e) 

• log_A_sum_init: calculate directly the logarithm of A_sum_init in case of overflow of the 
latter. 

• B_sum_init_PS-one: calculate - — — in the B(z) sum for the 1 — z linear transformation 

(1 - z)™ ^ ^ 

formula (see Eq. (l32i) ). 

• B_sum_init-PS -infinity: calculate in the B(z) sum for the - linear transformation 



formula (see Eq.(HT 



z 



cv-poly-der_tab_calc: calculate the coefficients of the derivative of the P{x) polynomial of 
SecdS]). 

cv-poly-der_calc: evaluate the derivative of the P{x) polynomial of Sec.(j6]) for a given x 
argument. 

min_n_calc: calculate the integer ric of Sec. ([6]) (A > 0) or return zero (A < 0) (see Sec. (EI) 
for method and notations). 

hyp_PS-zero: calculate F{z) for < 1 with Eqs.( l2]l3l) . 

hyp_PS-one: calculate F{z) for |1 — 2;| < 1 with the method of Sec.l lS.ip . 

hyp -PS -infinity: calculate F{z) for \z\ > 1 with the method of Sec. fl5.2p . 
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• hyp_PS_complex_plane_rest: calculate F{z) with a Taylor series expansion (see Eq. (H6|) ). 

• hyp_2Fl: calculate F{z) for arbitrary z with the method of Sec.©. 

• test_2Fl: test of F{z) accuracy from Eqs. fpi7|l . 

The routines meant to be used directly are hyp-2Fl and test-2Fl. hyp_2Fl demands as 
arguments the complex values of a, 6, c, and z in this order, and return the calculated 2F1 
function. To use test-2Fl, the tested complex value F{z) coming from hyp_2Fl has to follow a, 
6, c and z. A real number measuring its accuracy is returned, determined via Eqs. (l6|7|) . F\z) 
and F"{z) are therein calculated directly with hyp-2Fl. 

It is possible to use Gamma_inv-diff-eps, hyp_PS-zero, hyp_PS-one, and hyp -PS -infinity by 
themselves, but, for the three last functions, one has to pay attention to the modulus of their z- 
parameter. One must have < 1 and 1^1 > 1 respectively for hyp-PS-zero and hyp -PS -infinity. 
hyp-PS-one takes 1 — 2; as argument instead of z, so that the very close vicinity of z = 1 can 
be handled precisely. |1 — 2;] < 1 has to be verified as well. Moreover, it is necessary to ensure 
that 3fJ(c — a — 6) > in hyp-PS-one and 3f?(6 — a) > in hyp -PS -infinity. All other functions 
are not supposed to be handled by the user. 

Test examples are given in hyp -2F1 -example. cpp and hyp-2Fl-example.f90. It deals with a 
simple test of the 2F1 function, where the parameters a, 6, c, and z are given as inputs and the 
2F1 function evaluated value, as well the accuracy test of test-2Fl, are provided as outputs. 

10 The program PTG^wf 

The PTG potential and wave functions are programmed differently according to the C++ or 
Fortran 90 codes. In the C++ program, a class is defined for the potential, while a module is 
used for it in the Fortran 90 code. They are, however, very close so that the C++ routines will 
be delineated first, and one will just point out differences between both versions afterward. 

10.1 Routines of the C-\ — |- program 

All routines related to the PTG potential are situated in PTG-wf.cpp and PTG -wf -example, cpp 
The class defining the PTG potential PTG-class is built from a constructor demanding the 
following arguments in that order: 

• kin-fact: kinetic factor equal to , where mo is the rest mass of the particle. has 

the dimension of an energy times the square of a length. 

• /: integer angular momentum of the PTG potential. 

• Lambda, s,nu, a: parameters of the PTG potential defined in Sec. fl8.ip . 

This class possesses the main functions pertaining to the PTG potential formulas of Sec. fl8.ip : 

• Lambda2-sr-function: calculate the function A'^sr{y) given by Eqs.f H9|50i) . The y value 
must be given as argument. 
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y_search: calculate the y{r) function of Sec. (18.21) using Newton method and bisection. A 
value of one is returned if i/d > 0.99 (see Sec. (18. 21) ). 



• exp_minus-2x-iterative: calculate e with x = arctanh(?/) when y ~ 1 (see Sec. (l8.2p ). 
This function and the precedent take A^sr as argument. 

• effective-mass, effectivejmass-der: calculate the dimensionless effective mass /i(r) of the 
PTG potential (see Eg. (15 II) ). r is given as argument. Same for its derivative /i'(r). 

• V-cLcalc: calculate the ^-dependent Vf(r) potential part with Eqs. (l53ll56ll57f58p in energy- 
units. 

• operator(): calculate the PTG potential Vptg(^) of Eq.(!55l) in energy units. 

r always has the dimension of a length. Units can be chosen arbitrarily by the user and are 
invisible in the code. All arguments and returned values are double-precision real. Member 
functions Lamhda2_sr Junction and V-cLcalc are private as not intended to be used directly by 
the user. 

PTG wave functions and the associated test are calculated with non-member routines: 

• PTG-pole: calculate the wave function, derivative and second derivative of a parti- 
cle S-matrix pole state of principal quantum number n of a given PTG potential (see 
Eqs.( !79|80|) ). It can be bound, resonant, or antibound. Arguments are respectively: a 
reference on the PTG_class, the integer n, the number of points of the wave function, a 
table of real r radii, two complex numbers which will contain the Co and values of 
Eqs. (l90ll9T]) . and three tables of complex numbers where wave function and derivatives 
are stored. 

• PTG-Scat: Same for a scattering state of complex linear momentum k. Arguments are 
the same except for n which is replaced by fc, and a complex number is added between 

and wave function tables to store the value of C~ from Eq. fl9T|) . 

• PTGJ,est_calc: calculate a test of the PTG wave function accuracy from its Schrodinger 
equation (see Eq.lHHl)). Its arguments are: a reference on the PTG-class, the complex 
linear momentum k, a radius r and three complex numbers wielding the wave function, 
derivative and second derivative at r radius. 

Added to that, the complex function k_PTG-calc returns the complex linear momentum of a 
S-matrix pole state of principal quantum number n, taking as arguments the integers n and / 
and the real numbers Lambda, s, nu, a, and kinjact defining the PTG potential. Its unit is the 
inverse of a length. 

PTG_wf-example.cpp provides a test example of PTG wave functions. All parameters de- 
manded by the class PTG-class are input data, as well as the radial interval where PTG wave 
functions will be evaluated. The number of points considered in this interval are uniformly dis- 
tributed. Two states, a S-matrix pole and scattering state are calculated, so that the principal 
quantum number of the former and the linear momentum of the latter are demanded as inputs. 
The output consists of the evaluated PTG wave functions, first and second derivative, as well 
as in the PTG potential, effective mass and derivative at the same point r, and the accuracy 
test of PTGJ,est_calc. The input file has been configured so that the two possible formulas for 
calculating the PTG wave function Lp{r) in Eqs. (l79|80l) are utilized. 
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10.2 Fortran 90 program 

The Fortran 90 program routine names are the same, except for the operator () function of the 
class PTG-class which is renamed as V-PTG. They can be found in PTG-wf.fOO. Parameters of 
the PTG potential are initialized with the subroutine PTG-PARAMETERS, taking the same 
arguments as the constructor of the class PTG-dass. Arguments of other functions remain the 
same, except for the PTG_class member reference of C++ routines, absent in the Fortran 90 
code. 

11 Recommendations 

When utilizing hyp-2Fl, the main issue is the size of a, b, and c moduli. While the code provides 
virtually full numerical precision for all values of |a|, \b\, and |c| smaller or of the order of one, 
precision starts to deteriorate for larger moduli (see Tabs. (!3|4f5|) ). Spurious effects are the most 
important when the imaginary parts of parameters a, b, and c increase. However, the average 
error is smaller than the maximal error by several orders of magnitude for moderate and large 
\a\, \b\, and |c| (see Tabs. ( 151111131) ). This indicates that the instabilities are very localized in 
the a, b, c, and z four-dimensional complex space. Hence, calculations are usually meaningful 
for moderate |a|, \b\, and |c|. Testing functions is therein necessary to ensure the validity 
of the calculation, rendering, however, the time of calculation three times slower due to 2-^1 
calculations overhead. If accuracy is poor, it is possible to use directly hyp_PS-zero, hyp_PS-one, 
and hyp _PS -infinity along with transformations of Sec. ([3]) instead of hyp-2Fl, as mentioned in 
Sec. ([9]), as the method automatically chosen by the code to evaluate 2^1 may not be the most 
robust for some particular instances of a, b, c, and z. The polynomial case is almost devoid of 
numerical inaccuracy for polynomial degrees smaller than ten (see Tab. ([6])). 

If one calculates hyp-2Fl functions for integer-spaced a, 6, and c, the recurrence relations of 
Ref. [SU] could be helpful (see also Ref. [32] for a theoretical study of these relations). However, 
they become unstable if the calculated 2-F1 functions are minimal solutions of Eq.(jl]) [3H]- Hence, 
the calculation of 2-Pi(a + n,b + m,c + z) values, where (n, m, p) G Z, theoretically reachable 
from 2-Fi(a±l, 6, c; z), 2-P'i(o, 6±1, c; z) and 2-^1(0, 6, c±l; z) evaluations and recurrence relations 
[30] . would become too difficult to code. If the calculated 2-^1 function is minimal, backward 
recurrence with Miller's algorithm has to be used, while forward recurrence must be employed in 
the opposite case [SB] ■ As all three parameters a, 6, and c must be recurred, a general algorithm 
allowing the evaluation of 2Fi{a + n,b + m,c + p;z) would have to consider all possible forward 
and backward recurrences which can occur for arbitrary complex a, 6, and c, which would 
demand to consider too many particular cases. Moreover, they cannot solve the instability 
problem appearing for large |Q'[a]|, and |5J[c]|. They, may, however, prove useful for 

some particular problems. 

The situation concerning PTG wave functions is much better. Firstly, PTG pole states of 
principal quantum number n only involve 2-^1 polynomials of degree n (Jacobi polynomials) 
|24] . n rarely exceeds ten, so that the polynomial case of 2-F1 is therein sufficient. In case of the 
opposite, the Jacobi polynomials can be evaluated with stable recurrence relations. However, 
this situation is considered too particular to be included in our general code. Secondly, 2F1 
parameters of PTG scattering states will usually be of a small modulus, so that the aforemen- 
tioned instabilities will be absent in practice. 
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12 Conclusion 



The direct evaluation of the 2F1 function with power series and all its parameters and argu- 
ment complex has been attempted for the first time. The proposed method, relying on linear 
transformation theory and analytical properties of the Lanczos approximation of the Gamma 
function, has proved to be very efficient and allows to access all z values of the complex plane, 
the vicinity of negative integers b — a and c — a — b being rigorously treated. Its main drawback 
is its propensity to become unstable in general when |a|, |6|, and |c| increase. Such instabilities 
are, however, very localized, so that the code can still be used in most cases, as long as the 
accuracy of the functions has been verified to be satisfactory. Despite this obvious disadvan- 
tage, this code will surely allow to render possible calculations which before demanded the 
implementation of complicated algorithms for each encountered particular situation. 

The evaluation of PTG wave functions has been provided as a physical application, as 2-^1 
functions appear directly in their analytical formulas and because its implementation demands 
care due to its implicitly defined arguments. PTG wave functions have the advantage of being 
both analytical and bear bound, resonant and scattering states due to the finite range of 
the PTG potential. The published code permits a fast and stable evaluation of PTG wave 
functions, which will probably be of interest for studies in the domain of loosely bound and 
particle-emitting quantum systems. 
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Table 1: Relative accuracy of the 2F1 function obtained via Eqs. fl6|7|) for z = re* 3, with 
r = 0.99. The methods of Eqs. fl46ll47p . denoted as "ours", and the scheme of Ref. [33], denoted 
as "Biihring", are compared. \^{a,b,c)\ represents either |3fJ(a)|, \^{b)\ or |9ft(c)|, so that 0-1 
for \^{a,b,c)\ means that for x G {a,b,c}, < |9fJ(x)| < 1. For each line of the table, 30,000 
(a, b, c) sets of parameters are considered and are chosen randomly so that |9(a, b,c)\ < 1, with 
\Q{a,b,c)\ defined similarly to \^{a,b,c)\. Tmax is the maximal value obtained with Eqs. (!6H7I) 
and Tav is the average of all values given by Eqs. (l6|7|) . 



|3?(a,6,c)| 


Tmax (ours) 


Tmax (Biihring) 


Tav (ours) 


Tav (Biihring) 


0-1 


9.2 10-15 


1.4 10-11 


6.1 10-16 


4.7 10-15 


1-2 


9.7 10-15 


1.0 10-10 


1.1 10-15 


3.5 10-1^ 


2-5 


5.2 10-1^ 


5.1 10-10 


1.6 10-15 


1.6 10-13 


5-10 


9.5 10-13 


3.6 10-6 


6.7 10-15 


2.4 10-9 


10-15 


3.2 10-11 


1.1 10-2 


1.7 10-13 


2.9 10-5 



Table 2: Same as in Tab.([T]), but with r = 1.01. 



|3ft(a,6,c)| 


Tmax (ours) 


Tmax (Biihring) 


Tav (ours) 


Tav (Biihring) 


0-1 


9.0 10-13 


9.0 10-12 


3.0 10-15 


4.1 10-15 


1-2 


4.0 10-11 


9.2 10-11 


2.2 10-1^ 


3.6 10-1^ 


2-5 


9.6 10-11 


3.7 10-10 


6.8 10-1^ 


1.5 10-13 


5-10 


4.2 10-9 


3.0 10-6 


1.8 10-12 


1.5 10-9 


10-15 


8.5 10-^ 


6.9 10-2 


2.6 10-10 


1.6 10-5 
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Table 3: Relative accuracy of the 2-^1 function obtained via Eqs. fpiTl) . For each line of the table, 
10^ (a, b, c, z) sets of parameters are considered and are chosen randomly so that |^|oo < 3 and 
\'^{a, b,c)\ < 1. |3?(a, b,c)\, \'^{a, b,c)\, T^ax and Tav are defined similarly as in Tab.([T]). 



|3?(a,6, c)| 


T 

^ max 


T 


0-1 


1.0 10- 


12 


3.0 10- 


-16 


1-2 


4.1 10- 


11 


2.0 10- 


-15 


2-5 


3.1 10- 


10 


2.4 10- 


-14 


5-10 


5.0 10" 


-5 


8.8 10- 


-11 


10-15 


2.6 




1.2 10 


-5 



Table 4: Same as Tab. ([3]) except that the condition 1 < 10^(0,6,0)1 < 2 is verified. 



|3fJ(a,6, c)| 


T 

max 


T 

av 


0-1 


1.3 10- 


10 


7.7 10- 


-15 


1-2 


9.2 10- 


-9 


1.7 10- 


-13 


2-5 


3.3 10" 


-5 


6.6 10- 


-10 


5-10 


8.7 10" 


-3 


2.2 10" 


-8 


10-15 


5.8 




4.4 10' 


-5 



Table 5: Same as Tab. ([3]) except that the condition 2 < \Q{a,b,c)\ < 5 is verified. 



\^{a,b,c)\ 


T 

max 


T 

av 


0-1 


7.0 10-2 


3.4 10-^ 


1-2 


1.5 10-1 


2.7 10-^ 


2-5 


7.3 


1.5 10-^ 
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Table 6: Relative accuracy of 2-^1 polynomials obtained by way of Eqs. ipiTI) . 2^1 functions 
are polynomials of degree running from zero to ten, so that a is fixed to minus the degree. 
Parametrization is the same as in Tab. ([3]), except that 3fJ(6, c) and c) are used instead of 
3fJ(a, 6, c) and $5(0,6, c). One has here c)| < 10. 



mb,c)\ 


T 

max 


T 


0-1 


5.1 10-*^ 


1.2 10^ 


-12 


1-2 


8.1 10-^ 


2.3 10- 


-12 


2-5 


4.1 10-*^ 


2.3 10- 


-12 


5-10 


3.7 10-^ 


8.9 10- 


-12 


10-15 


3.8 10-6 


3.9 10- 


-11 
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